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The thermodynamical properties are calculated for a three-dimensional model of N harmonically 
interacting spin-polarized fermions in a parabolic potential well. The obtained dependences of the 
chemical potential and of the internal energy on the complete range of the temperature and of 
the number of particles turn out to obey a scaling law, similar to the scaling from the continuum 
approximation for the density of states. The calculation technique is based on our path integral 
approach of symmetrized density matrices for identical particles in a parabolic confining well. 
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I. INTRODUCTION 

In two preceding papers the present authors extended the method of symmetrized density matrices to systems 
confined in a parabolic well |jj (hereafter referred to as I) and used this method to obtain expressions for the density 
and the pair correlation function Q (hereafter referred to as II). The evaluation of the internal energy, the specific 
heat, the moment of inertia fipj, the density and the static response functions was performed for bosons, inspired by 
the recently observed Bose-Einstein condensation and the theoretical work around this phenomenon using other 
methods f| 0. In both papers I and II the general expressions for most of the quantities mentioned above are also 
given for Fermi-Dirac statistics. In the present paper a method is presented to explicitly evaluate the thermodynamic 
quantities for spin-polarized fermions. The model is a parabolic well containing N fermions all in the same spin state 
and interacting through a harmonic two-body potential that may be either attractive or repulsive. 

A quantum dot would be a physical system that could be described by such a model, if also a magnetic field 
is taken into account to freeze away the opposite spin states. When no magnetic field is taken into account two 
spin states (spin up and spin down) should be present in the model. Recently Jll| it was proposed to investigate 
also confined fermions in the same experimental configuration used for the bose-alkali metals. The Thomas-Fermi 
approximation |l6| ] was used to study the spatial distribution of these trapped fermion systems. In order to analyze 
these more physically relevant models in the future, we first develop in the present communication the required basic 
techniques for spin-polarized fermions. The model has also some importance in itself, because it can be used to test 
new approaches to Monte Carlo simulations of interacting fermions such as many-body diffusion |T^-|l9[]. 

The paper is organized as follows. In section II, we collect the expressions from I and II for the fermion case, and 
we show in section III how the chemical potential, the free energy and the internal energy can be obtained for a given 
number of fermions as a function of the temperature. Subsequently the low temperature limit is considered, and the 
groundstate energy is evaluated in section IV. In the last section a discussion and the conclusions are given. 

II. FERMION OSCILLATORS 

In this section the basic formulas which have been derived in the path-integral treatment of I and II are summarized 
and rewritten in such a way that they are more appropriate for dealing with fermions, in particular for the numerical 
treatment. Before doing so, it is instructive to point out where the numerical accuracy problems are coming from. 
Having pinpointed their origin, a method is proposed to accurately evaluate the relevant thermodynamic quantities. 
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A. Summary of previous results 



We consider a model of N fermions with parallel spin in an harmonic confinement potential, and with a quadratic 
interparticle interaction. The one-body potential energy V± and the two-body potential energy V2 of the model system 
are given by 

02^ 2 N 

j=i 3,1=1 

The two-body interaction is assumed to be repulsive; replacing — lo 2 by lo 2 in V2 gives the case of attraction. In each 
dimension we found one degree of freedom (the center of mass) with frequency tt, and TV — 1 degrees of freedom with 
frequency w given by 



w 



Vtt 2 - Nuj 2 , (2.2) 



which means that the frequency fl of the center of mass is larger than the frequencies w of the degrees of freedom in 
the relative coordinate system. Changing the sign of uj 2 allows to obtain the case with w larger than f2. 

In our path-integral treatment presented in I, a recurrence relation was obtained for the partition function Z/ (N) 
corresponding to the degrees of freedom with frequency w in the relative-coordinate system. Introducing 

b = e-P hw (2.3) 

for brevity in the notations, we found that: 

m=0 

This recurrence relation applies for bosons (£ = +1) and for fermions (£ = —1). The subscript / refers to identical 
particles, which can be specified to be fermions (subscript F) or bosons (subscript B). The total partition function 
Zj (N) only differs from Zj (N) by a factor which accounts for the center-of-mass contribution 

. . / sinh i Bw \ 3 , , , 

*W-(ssfei) z ' (Ar) - (2 - 5) 




B. The "sign"— problem and the canonical ensemble 



For three-dimensional (3D) fermions, the contribution (2.4) to the partition function from the relative degrees of 



freedom clearly illustrates the kind of numerical inaccuracies which originate for the fermion case £ = — 1. If the 
partition functions for 1, 2, (AT — I) particles are known, Cramer's rule can be used to calculate the partition 



function for N fermions. Factorizing the denominators in the partition function (2.4) by introducing the quantities 

JV 

z N = b-i N Z F (N)~[[(l-V) 3 , (2.6) 

i=i 

a careful analysis shows that zn are polynomials in b. Typical terms of the expansion in powers of b are summarized 
in table |. In table || the polynomial zio is given in full detail. 



These expressions clearly illustrate that the recurrence relation (|2.4| ) with its alternating signs is numerically not 
stable: the leading terms are of order b M where M increases drastically if the number of fermions increases. Never- 
theless, the expression for zio, e.g., is useful to check expressions for the partition function or derived quantities on 
their accuracy. Because for the fermions, unlike the boson case, solving the recurrence relations thus runs into severe 
numerical problems, we will use the generating function technique for the actual calculation of the free energy and 
the internal energy. To convince ourselves that numerical inaccuracies have been avoided, the internal energy of the 
model for up to 10 fermions has been calculated both ways, i.e. from the recurrence relations and with the generating 
function technique, and the results of both methods coincide. How the chemical potential and the internal energy are 
calculated will be elaborated in the next section. 
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III. THE THERMODYNAMIC PROPERTIES 



The generating function 3^ (u) corresponding to the partition functions 7Lp (N) is defined in the standard way as 

oo 

(it) = ^ u N Z F {N) . (3.1) 



-F ' 



JV=0 



As shown in I, it is given by 

E F (u) = exp -}'- ) 1 ■ (3.2i 




This means that in our model the internal degrees of freedom are represented by a system of non-interacting oscillators 
with frequency w. 5f (u) is then formally the grand-canonical partition function of that subsystem. However, it is 
not the grand-canonical partition function of the full model system with interaction for two reasons: first one has to 
take the center-of-mass correction into account, and second the eigenfrequency w in the relative coordinate system 
depends on the number of particles. But given w the full mechanism of the generating functions is applicable in the 
relative coordinate system, provided afterwards the necessary center-of-mass corrections are taken into account. 

The partition function Zf (N) from the internal degrees of freedom can be obtained by inverting the defining Taylor 
series ( pTl] ) 

z *w = hliF&* t (3 - 3) 

where C is a closed contour in the complex z plane around the origin. The generating function 5_p (z) is unaccessible 
for numerical purposes. However, considering a circular contour with radius u one obtains 



Z F (N) = — \^ Ef y ) e- im d8 = — /^exp [ln3 F (ue ie ) -Nlnu] e~ iNe dO. 



(3.4) 



The extremum of [lnSjr (ue l6 ) — Nlnu] on the real axis satisfies the condition N = u-^ LnSj? (u) . Using ( |3.2D this 
requirement becomes 

v— ^ 

which is precisely the result which one would obtain from the "grand canonical" treatment with u = e' 3 ^, taking into 
account the degeneracy | {v + 1) {v + 2) of the ^th energy level. Factorizing out the steepest descent contribution 
'E.f (u) /u N obtained this way, one finds 

Z F (N) ==^P- r*(0)d0, (3.6) 
u Jo 

9(0)=- 3 i^ a ) e- iN9 , (3.7) 
7T aj? (u) 

where (8) is a real function, suitable for numerical integration if u = is determined. 

The advantage of a procedure based on the generating function is that all contributions to 5 F (u) turn out to be 



positive, in contrast to the direct determination of the partition function (2.4) which numerically involves severe sign 
problems, as argued in the previous section. 



A. The chemical potential 



The chemical potential has to be determined from the requirement (3.5). There are clearly two cases to be consid- 



ered. For sufficiently low temperature, \x will be larger than ~Hw, but at high temperature /i might be smaller than 
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For the case fi > § Hw, the behavior of the denominator is fundamentally different for the energy levels t v < \i as 
compared to those satisfying e u > /i. Assume that fi will be found in the interval between two consecutive 

levels 

/Lt = cl — hwa; a € [0, 1[, (3-8) 

and treat the levels e„ < /i separately from those with e v > /i. We start the discussion under the assumption L > 1 
and rewrite AT as: 

A 1 (L-i/+1)(L-i/ + 2) 1 (L + 1) (£ + 2) y> 1 (L + i/ + 1) (L + t/ + 2) b a+v 

- 2 1 + 6" a+l/ + 2 1 + 6 a + ^ 2 1 + b a+lJ 

i/=l w=i 

Let A^ denote the (not necessarily integer) value which the right hand side of this equation would take for a = 0, i.e. 
if the chemical potential would exactly be located at e.^. One finds after some algebra (provided L > 0) 

N l = 1 -(l + 3 -)(L + 2)(L + 1) + (2L + 3)± J"L+ V 1(^ + ^1)^ + ^2)6-1 
6\ 2/ V M ; v ^1 + 6" ^2 1 + 6^ 6 

v 7 i/=l y=L+l 

1/3 

The result L < (6Nl) is important for developing a sure-fire algorithm to find the chemical potential. If one starts 

1 /3 

from a value L which is the largest integer smaller than (6N) 1 , it implies that /i < e^. One can then decrease the 
value of L step by step, until a value of L is found for which Nl > N but Nl-i < N. 

Fortunately, this procedure can also be used if /i is smaller than t he lowest energy of the model. It then results in 
a negative value of L, which is related to the chemical potential in ( |3.8[ ) by formally filling out ez, with L negative. 

A fresh refining routine is then started to find a € [0, 1[, which is bracketed as required for sure-fire root finding 
programs. The actual determination of a (for both cases /i > ^hw and fi < ^hw) is straightforward from the equation 
for A~ if written in the appropriate form for numerical treatment: 

^2 1 + b L - v - a ^ 2 v M ; I , b L+v+a ^ else v ' 

The determination of the chemical potential along these lines presents no numerical difficulties. Because the Fermi 
energy is of order hw (QN) 1 ^ 3 , it seems natural to express kT in units of hwN 1 ^ . This scaling factor turns out to 
be surprisingly good, as is shown in Fig. 1 where the temperature dependence of the chemical potential /i (T) in 
units of /Lt (0) is plotted against t = ft j^i/s - Introducing a density of states and making the continuum approximation 
(ksT -C Hw) the corresponding scaling law |l5| , [l6|] with the Fermi energy is implicit. 



B. The free energy 

Having determined the chemical potential, the free energy Fp (N) — — i lnZ^ (AT) can be evaluated from 

¥ F (N) = 4 0) (N) - i In ( [ *(0)d9), (3.10) 



^ } (N)=-~ln^l, (3.11) 



where (N) is the zero-order steepest descent result, which would be obtained from the "grand-canonical" treat- 
ment. The correction factor involving \& (0) accounts for the finite number of particles. For reference, 'J (9) is shown 
for AT = 10 in Fig. 2 as a function of 9 for various values of the temperature. If A?" increases, 'J (9) becomes increasingly 
concentrated near the origin 9 = for nonzero temperatures. 

The resulting free energy per particle as a function of temperature is shown in Fig. 3-5, for N = 1, N = 10 and 

1 /3 

A" = 100 respectively, in units of Hw (6N) ' proportional to the Fermi energy. For comparison, the contribution 

Fp?^ (N) is also plotted (dashed lines) . As expected, this steepest descent contribution becomes increasingly accurate 
if the number of fermions increases. 



4 



C. The internal energy 



The contribution of the relative degrees of freedom to the internal energy 

Uf = lp {(3¥f) = ¥f ~ T ^ ¥f (3 ' 12) 
can be obtained from the free energy obtained above by numerical differentiation. The internal energy per particle is 

1 /3 

plotted in Fig. 6 in units of hw (67V) ' proportional to the Fermi energy. A similar scaling law is observed as for the 
chemical potential. 

In one dimension the thermodynamical properties of harmonically interacting bosons and fermions can be derived 
from each other. This case has been studied in I. In higher dimensions, the fermion internal energy is smoothly 
decreasing with decreasing temperature, whereas for the boson case it shows sudden changes in slope, related to the 
condensation. 

To within the numerical accuracy, the results are in agreement with the standard description from the generating 
function treatment using 

^ ^l („ + l)(„ + 2) fr + 3/2) 

^e„n v = hw^- -——j . (3.13) 

These results can also be compared to the internal energy Ufr rcc (TV) which one would obtain from the recurrence 
relation. In terms of the expressions for discussed above one then obtains 

This calculation is in practice only feasible for a limited number of particles TV < 10, and for these cases it coincides 
within the numerical accuracy with the results plotted in Fig. 6. 



IV. THE GROUNDSTATE ENERGY 



In this section the low-temperature limit will be considered. By counting the number of occupied energy levels, 
the dominant contribution to the partition function can then easily be calculated, taking into account the degeneracy 
\ (y + 1) (v + 2) of the levels with energy e v — hw (y + I). The calculation is done first with the Fermi level L fully 
occupied. The number of particles Np required for this assumption to hold is 

L 1 1 

Nf = E 2 {V + X) {V + 2) = 6 {L + 3) {L + 2) (L + 1] ■ (41) 

v— 

1/3 

Consequently, the Fermi energy is of order (67V) ' hw. The total energy Uf associated to the case with the L-th level 
fully occupied is 

U. 



^1^ + 1)^ + 2)^^+^ =hw(L + 3)(L + l)(L + 2f. (4.2) 



v=0 



For a limited number of particles, the number of particles and the energy Up* are shown in Table III. 

For an arbitrary number TV of fermions not filling the Fermi level completely, the determination of the ground state 
ener gy is slightly more involved. We first determine the number of particles Np < N which fills the level L. From 



(hi) it follows that the highest fully occupied level L is given by 

1 . \-V3 \ 

(4.3) 



// 1 \V3 ,/ -, X -1/3 

L = integer 3iV + - \/3 6 N 2 - 3 + - 37V + - a/3 6 7V 2 - 3 



9 v J 3 V 9 

From this level L, one can determine the corresponding TVp and Uf 



5 



The resulting formula for the leading term in the partition function for N particles in the zero temperature limit is 
b 3N/2 b M N j which defines the power M N as 

M N = (L + 1) (n - i (L + 2) (L + 3) (L + 4)) . (4.4) 

The ground state energy Eq with N particles is Eq = Up + (N — Np) Hw (L + 1 + |) because the remainder of the 
particles is in the level L + l, and consequently 

E = (L + 4) (L + 3) (L + 2) (L + 1) + Nhw ( L + 1 + | J . (4.5) 



V. CONCLUSION AND DISCUSSION 



In this communication we have given a short review of the calculation techniques for fermions, described in I and 
II for identical particles in general. Next a numerical analysis of the chemical potential and of the free energy is 
made for a given expectation value of the number of particles as a function of temperature. In this analysis, we could 
easily illustrate what the consequences are of the minus sign coming from the anti-symmetric representation of the 
permutation group in the expression for the partition function. Even when the expressions are known analytically, the 
plot of a relatively smooth function such as the free energy requires special techniques as a consequence of numerical 
instabilities due to a sign problem in the recurrence relation for the partition function. The necessity of such techniques 
can be checked by attempting a calculation of a few limits, which lead to the application of De l'Hopital's rule many 
times, even proportional to the square of the number of particles in the system. 

It should be noted that the model contains only spin-polarized fermions. In quantum dots, its production would 
require a magnetic field. In our model this field is not included. However as we have shown in I, the expressions 
for the partition function for fermions in the presence of an external magnetic field can be obtained with the same 
calculation technique. The influence of the magnetic field on the chemical potential and the specific heat of our model 
has not been studied yet for fermions. In alkali metal vapors the spin polarization of the fermions would be inherent 
to the experimental technique |20| . 

The chemical potential and the internal energy exhibit a scaling law in the sense that it is an almost universal 
function of the temperature when plotted in the indicated scaled units. Although we strongly suspect that the scaling 
comes via the Fermi level of the confined system as is the case in the continuum limit, we have no mathematical proof 
of this observation in the low temperature case fegT -C hw. 

We did not compare the present approach with other theories using the same or an analogous model. It should 
however be stressed that as far as we could check we presented here the first results obtained with a new scheme for 
the evaluation of path integrals for fermions. 
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Z 2 = 6(3+6 2 ) 

z:i = b 2 (3 + 106 + 6b 2 + 6b 3 + 7b 4 + 36 s + b r ) 

24 = b 3 (1 + 156 + 276 2 + 626 3 + 636 4 + 876 5 + 806 6 + 876 7 + O (6 8 )) 

25 = 6 5 (6 + 376 + 1056 2 + 2316 3 + 4136 4 + 6696 s + 9216 6 + 11976 7 + O (6 8 )) 

z 6 = 6 7 (15 + 756 + 2906 2 + 6876 3 + 15906 4 + 29946 s + 53046 6 + 83886 7 + O (6 8 )) 

27 = 6 9 (20 + 1356 + 5436 2 + 16456 3 + 42 066 4 + 938 1 6 s + 19131b 6 + 358026 7 + O (6 8 )) 

2 8 = 6 11 (15 + 1736 + 7806 2 + 28716 3 + 82 966 4 + 2 14536 s + 491106 6 + 1047236 7 + O (6 8 )) 

29 = 6 13 (6 + 1356 + 84 76 2 + 36 1 26 3 + 123 486 4 + 36 1 666 s + 939726 s + 2235726 7 + O (6 8 )) 
210 = 6 1S (1 + 576 + 6 1 56 2 + 326 1 6 3 + 135 036 4 + 453 456 s + 1346106 6 + 3579336 7 + O (ft 8 )) 

TABLE I. Reduced partition function z N for A = 2, 3, 10. 



2io/6 ls = 1 + 576 + 6156 2 + 326 1 6 3 + 135 036 4 + 453 456 s + 1346106 6 + 3579336 7 + 8790546 s + 20106846 9 + 43451286 10 
+89180286 11 + 175221216 12 + 330747666 13 + 602694756 14 + 1062918456 15 + 1820052216 16 + 3031594506 17 
+4922982736 18 + 7805097696 19 + 12101 169696 20 + 18367 968086 21 + 27328 288896 22 + 39890 231586 23 
+57179 095546 24 + 80544 274896 2s + 1 11580 118886 26 + 1 52 106 1 58 466 27 + 2 04 163 941636 28 + 2 6 9 9 5 5 4 65006 29 
+3 51805 186 266 30 + 4 520 4 5 9 1 7716 31 + 5 729 4 3 2 73366 32 + 7 165 2 9 671646 33 + 8 8 4 50 1 887056 34 + 10 780 1 9 8 3926 1 
+12 97604 687676 36 + 15 42935 025606 37 + 18 12786 766656 38 + 21 04869 193096 39 + 24 15851 836596 40 
+2741281 606566 41 + 30 75731 389756 42 + 34 12822 384386 43 + 3745505 135706 44 + 40 66182 875296 4s 
+43 67101 708776 46 + 46 40538 232416 47 + 48 79244 290706 48 + 50 76635 376066 49 + 52 27216 273326 50 
+53 26698 445206 s1 + 53 72332 781856 s2 + 53 62897 297986 s3 + 52 98885 994926 s4 + 51 82338 664596 s5 
+50 16872 645286 s6 + 48 07367 472156 s7 + 45 59865 530106 s8 + 42 81161 407776 s9 + 39 786 2 8 5 8326 60 
+36 59783 309826 61 + 33 32107 425996 62 + 30 02657 983866 63 + 26 77945 028256 64 + 23 63641 326726 65 
+20 6 4551 180 986 66 + 17 8 4 4 4 2 4 8 2 226 67 + 15 26113 4 00 1 96 68 + 12 9 133 9 2 1 7396 69 + 10 8 1 16 3 79486 70 
+8 95194 774036 71 + 7 33263 427906 72 + 5 94028 454766 73 + 4 75902 667656 74 + 3 76992 180576 75 + 2 95259 676696 7f 
+2 28593 937966 77 + 1 74930 025806 78 + 1 32289 338576 79 + 98853 842676 80 + 72976 285726 81 + 53215 147266 82 
+38322 461676 83 + 27251 033406 84 + 19129 722246 8S + 132 5 5 79226 86 + 90 6 2 9 70266 87 + 6114143556 88 
+4068374556 89 + 2669987846 90 + 1727465286 91 + 1101898176 92 + 692571696 93 + 429000696 94 + 261701526 95 
+157286646 96 + 93038946 97 + 54211806 98 + 31067466 99 + 17537556 100 + 9728196 101 + 5317556 102 
+2852046 103 + 1509036 104 + 781976 105 + 400386 106 + 200106 107 + 99286 108 + 47586 109 + 22986 110 
+10566 111 + 4926 112 + 2136 113 + 996 114 + 396 lls + 186 116 + 76 117 + 36 118 + 6 120 

TABLE II. Reduced partition function zio- 
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VF 


iV F 


Uf /fow 


VF 


N F 


Uf /hw 


VF 


N F 


Uf /hw 


1 


4 


9 


10 


286 


2574 


60 


39711 


3693123 
2 


2 


10 


30 


15 


816 


10404 


70 


62196 


3358584 


3 


20 


75 


20 


1771 


58443 
2 


80 


91881 


11301363 
2 


4 


35 


315 
2 


25 


3276 


66339 


90 


129766 


8953854 


5 


56 


294 


30 


5456 


130944 


100 


176851 


27058203 
2 


6 


84 


504 


35 


8436 


234099 


150 


585276 


66721464 


7 


120 


810 


40 


12341 


777483 


200 


1373701 


416231403 


8 


165 


2475 


45 


17296 


609684 


250 


2667126 


504086814 


9 220 


1815 


50 


23426 


913614 


300 


4590551 


20795 19603 
2 


10 


286 


2574 


55 


30856 


1319094 


400 


10827401 


65289 22803 
2 



TABLE III. Number of particles Nf and total energy Uf with the energy level vf fully occupied. 

Figure captions 

Fig. 1: Scaled chemical potential jffjQx as a function of the scaled temperature t = h ^!/ 3 for 10, 100, 1000 and 
10000 fermions. For reference, this quantity is also plotted for 1 particle. 

Fig. 2: The integrand '5 (9) of equation ( |3.7j ) for 10 fermions as a function of 9 for various values of the temperature, 
expressed in units of T rc f = hwN 1 / 3, /k. 

Fig. 3: Scaled free energy per particle / = h ^^i/ 3 as a function of the scaled temperature t = h ^ 1/3 for 1 particle. 
For comparison, the zero-order steepest descent contribution is also plotted (dashed line). 

Fig. 4: Same as Fig. 3 but for 10 fermions. 

Fig. 5: Same as Fig. 3 and 4 but for 100 fermions. 

Fig. 6: Scaled internal energy per particle u = h ^^i/3 as a function of the scaled temperature t = h ^ 1/3 for 10, 
100, 1000 and 10000 fermions. For reference, the result for 1 particle is also plotted. 
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